

clear all
clc

b = str2double(getenv('SLURM_ARRAY_TASK_ID'));

load('Match_EndShare_minmax_small.mat')

cf = 0;

datsize = 100;
S       = 100;
B       = 100;
out     = 1;

%% Specify Model Characteristics

ModelCt = 19;

maxJJ_NN = 21;
maxJJNN = 98;
maxNN = 15;
maxN = size(data,1);
maxM = 426;
maxN0 = size(data0,1);
J = 27;

%% Counterfactual

royal25     = 0;
if cf == 2 || cf == 3 || cf == 6 || cf == 7 || cf == 8 || cf == 9 || cf == 10 
royal25     = 1;
end

%% Halton draws (use a different set of S draws for each MC)

StateVec = randi(1e9,[B,2,S]);
v = b;

FirmError_Sim = zeros(maxN,S);
ParcelError_Sim = zeros(maxN,S);

for s = 1:S

const.NObs = maxN; 

rng(StateVec(v,1,s))
skip = randi(1e3);
leap = randi(1e2);

p = haltonset(1,'Skip',skip,'Leap',leap);
p = scramble(p,'RR2');
draw = net(p,const.NObs);
draw = norminv(draw);

% FirmError_Sim(:,s) = draw(firmID,:);
FirmError_Sim(:,s) = draw;

clear draw p

const.NObs = maxN; 

rng(StateVec(v,2,s))
skip = randi(1e3);
leap = randi(1e2);

p = haltonset(1,'Skip',skip,'Leap',leap);
p = scramble(p,'RR2');
draw = net(p,const.NObs);
draw = norminv(draw);

ParcelError_Sim(:,s) = draw;

clear draw p

end

%% Equilibrium Convergence

SampleMMean_mod     = zeros(maxM,maxJJ_NN,ModelCt + 1);
SampleMVar_mod      = zeros(maxM,maxJJ_NN,ModelCt + 1);
SampleMCov_mod      = zeros(maxM,maxJJNN,ModelCt + 1);

SampleJMean_mod     = zeros(maxM,maxNN,ModelCt + 1);
SampleJCov_mod      = zeros(maxM,maxJJNN,ModelCt + 1);
SampleJVar_mod      = zeros(maxM,maxNN,ModelCt + 1);
SampleJDist_mod     = zeros(maxM,maxJJNN,ModelCt + 1);

SampleJntDist_mod   = zeros(maxJJNN,ModelCt + 1);
SampleCov_mod       = zeros(maxJJNN,ModelCt + 1);

Share_mod           = zeros(maxM*J,ModelCt + 1);
ShareMean_mod       = zeros(maxM,maxJJ_NN,ModelCt + 1);
ShareFirm_mod       = zeros(maxM,ModelCt + 1);
ShareCov_mod        = zeros(maxM,maxJJ_NN,ModelCt + 1);

shareOUTCF_st_mod   = zeros(maxN,ModelCt + 1);
QualityCF_mod       = zeros(maxN,ModelCt + 1);

ParcelUtil_mod      = zeros(ModelCt+1,maxN0);
FirmUtil_mod        = zeros(ModelCt+1,maxN0);
Quality_mod         = zeros(ModelCt+1,maxN0);
ParcelNotMatched_mod = zeros(maxN0,5,ModelCt+1);

landsqft_mod        = zeros(ModelCt+1,maxN0);
landsqft_notM_mod   = zeros(ModelCt+1,maxN0);
gas_mod             = zeros(ModelCt+1,maxN0);
gas_notM_mod        = zeros(ModelCt+1,maxN0);

LandByMarket = zeros(maxM,2,S,ModelCt+1);

rankA_valns_mod     = zeros(J,maxN0,ModelCt+1);
qualA_valns_mod     = zeros(J,maxN0,ModelCt+1);
rankA_val_mod       = zeros(J,maxN0,ModelCt+1);
qualA_val_mod       = zeros(J,maxN0,ModelCt+1);

diffnsMean          = zeros(maxN-1,ModelCt+1);
diffMean            = zeros(maxN-1,ModelCt+1);

for v = 0:ModelCt
%cd C:\Users\Ashley Vissing\Desktop\matlabhold\SLURM_match2018\ 
    check = exist(['DataResult_v' num2str(v) '_ival' num2str(b) '_myopic_out_' num2str(datsize) '.mat'],'file');
    if check == 2
        load(['DataResult_v' num2str(v) '_ival' num2str(b) '_myopic_out_' num2str(datsize) '.mat'], 'betaStar' )
        load(['DataResult_v' num2str(v) '_ival' num2str(b) '_myopic_out_' num2str(datsize) '.mat'], 'theta_score')
        load(['DataResult_v' num2str(v) '_ival' num2str(b) '_myopic_out_' num2str(datsize) '.mat'], 'hat' )
        load(['DataResult_v' num2str(v) '_ival' num2str(b) '_myopic_out_' num2str(datsize) '.mat'], 'bonus_hat')
        
        load('Match_EndShare_minmax_small.mat')

        %% Define data matricies and parameter counts

        N = size(data,1);
        propID = data(:,1);
        firmID = data(:,5);

        % Data needs to be sorted by market
        market = data0(:,2);
        marketindex = unique(market);
        marketindex = [(1:size(marketindex,1))' marketindex];

        M = max(marketindex(:,1));
        MarketSizeVec = zeros(M,1);
        for m = 1:M
            row = find(market == marketindex(m,2));
            MarketSizeVec(m,1) = size(row,1);
        end
        J = unique(data(:,5)); J = size(J,1);

        FirmMatch = data(:,5);
        FirmMatch0 = data0(:,5);

        HOLD = [data0(:,5) data0(:,8)];
        FirmMatchConstraint = [(1:J)' zeros(J,1)];
        for j = 1:J
           row = find(HOLD(:,1) == j);
           FirmMatchConstraint(j,2) = size(row,1);
        end
        clear HOLD
        FirmMatchConstraint = FirmMatchConstraint(:,2);

        %% Share data

        % Stacked share of leases signed by a firm in a wellpad
        share0 = data(:,7)./data(:,97);
        share0(isnan(share0)) = 0;

        % Firm share observed in the data to match in the algorithm
        FirmShare0 = [data(:,2) data(:,5) data(:,7) data(:,97)];
        FirmShare0 = unique(FirmShare0,'rows');
        FirmShare0 = FirmShare0(:,3)./FirmShare0(:,4);
        FirmShare0(isnan(FirmShare0)) = 0;
        
        %% Data vectors

        pipeDist        = data(:,49); pipeDist0        = data0(:,49);
        wellDist        = data(:,50); wellDist0        = data0(:,50);
        landSize        = data(:,48); landSize0        = data0(:,48);
        prodwell        = data(:,27);  prodwell0       = data0(:,27);
        enforce         = data(:,28);  enforce0        = data0(:,28);
        app_val         = data(:,17);  app_val0        = data0(:,17);
        parcelToWP      = data(:,57);  parcelToWP0     = data0(:,57);
        parcelToApi     = data(:,60);  parcelToApi0    = data0(:,60);
        parcelInt       = data(:,16);  parcelInt0      = data0(:,16);
        landman         = data(:,33);  landman0        = data0(:,33);
        largeOp         = data(:,21);  largeOp0        = data0(:,21);

        hat_ind = 0;

        % Imputted bonuses
        if bonus_hat == 1
            bonus_ind = 1;
        else
            bonus_ind = 0;
        end

        % Lease quality
        Clause  = data(:,67 + hat_ind); Clause0 = data0(:,67 + hat_ind); 

        % Well profit
        expFutureParcel = data(:,53 + bonus_ind);  expFutureParcel0  = data0(:,53 + bonus_ind);
        expFutureFirm   = data(:,51 + bonus_ind);  expFutureFirm0  = data0(:,51 + bonus_ind);


        %% Adjust the royalty
        if cf > 0
        
        %% Add clauses

        if cf == 3
            % Full set of clauses
            Clause  = repmat(max(Clause),[N,1]);
            Clause0 = repmat(max(Clause),[N/J,1]);
        %     Clause  = repmat(1,[N,1]);
        %     Clause0 = repmat(1,[N/J,1]);
        end
        if cf == 6
           % Add a clause
           Clause = data(:,67);  Clause0 = data0(:,67);
        end
        if cf == 7
           % Add env clause
           Clause = data(:,68);  Clause0 = data0(:,68);
        end
        if cf == 8
           % Add noise clause
           Clause = data(:,69);  Clause0 = data0(:,69);
        end
        if cf == 9
           % Add subsurface easement clause
           Clause = data(:,70);  Clause0 = data0(:,70);
        end
        if cf == 10
           % Add surface damage clause
           Clause = data(:,71);  Clause0 = data0(:,71);
        end

        ParcelIndex     = [expFutureFirm parcelToWP landSize app_val pipeDist wellDist parcelInt Clause];
        Operator        = [landSize.*largeOp wellDist.*largeOp parcelInt.*largeOp];
        Landman         = [landSize.*landman wellDist.*landman parcelInt.*landman];

        ParcelIndex     = [ParcelIndex Landman Operator];

        ParcelIndex0    = [expFutureFirm0 parcelToWP0 landSize0 app_val0 pipeDist0 wellDist0 parcelInt0 Clause0];
        Operator0       = [landSize0.*largeOp0 wellDist0.*largeOp0 parcelInt0.*largeOp0];
        Landman0        = [landSize0.*landman0 wellDist0.*landman0 parcelInt0.*landman0];

        ParcelIndex0    = [ParcelIndex0 Landman0 Operator0];

        FirmMatrix      = [expFutureParcel  prodwell  enforce  wellDist.*largeOp   landSize.*largeOp  Clause];
        FirmMatrix0     = [expFutureParcel0 prodwell0 enforce0 wellDist0.*largeOp0 landSize0.*largeOp0 Clause0];


        NN = size(ParcelIndex,2);
        JJ = size(FirmMatrix,2);

        % Parcel preferences
        gamma = [1 betaStar(1,1:NN - 1)];
        delta = [1 betaStar(1,NN:NN + JJ - 2)];
        alpha = betaStar(1,NN + JJ - 1);

        % Firm Value of parcels
        ParcelUtility = ParcelIndex*gamma';
        FirmUtility = FirmMatrix*delta';

        %% Counter-factual
        
        SampleMMean_cf       = zeros(M,JJ + NN,S);
        SampleMVar_cf        = zeros(M,JJ + NN,S);
        SampleMCov_cf        = zeros(M,JJ*NN,S);

        SampleJMean_cf       = zeros(J,NN,S);
        SampleJCov_cf        = zeros(J,NN*JJ,S);
        SampleJVar_cf        = zeros(J,NN,S);
        SampleJDist_cf       = zeros(J,NN*JJ,S);

        SampleJntDist_cf     = zeros(NN*JJ,S);
        SampleCov_cf         = zeros(NN*JJ,S);

        Share_cf             = zeros(M*J,S);
        ShareMean_cf         = zeros(M,NN + JJ,S);
        ShareFirm_cf         = zeros(M,S);
        ShareCov_cf          = zeros(M,NN + JJ,S);

        FirmMatchCF   = zeros(S,1);
        shareOUTCF_st = zeros(N,S);
        QualityCF_st  = zeros(N,S);

        ParcelUtil_sim  = zeros(S,N/J);
        FirmUtil_sim    = zeros(S,N/J);
        Quality_sim     = zeros(S,N/J);
        ParcelNotMatched_sim = zeros(S,N/J,5);

        landsqft = reshape(landSize,[J,N/J]);
        gas      = reshape(data(:,49),[J,N/J]);

        landsqft_sim    = zeros(S,N/J); landsqft_notM_sim = zeros(S,N/J);
        gas_sim         = zeros(S,N/J); gas_notM_sim = zeros(S,N/J);
        Quality         = Clause;

        % Market indexes - stacked
        marketUni = unique(market);
        MarketMat = [];
        for m = 1:M
            row = find(market == marketUni(m,1));
            MarketMat = [MarketMat ones(1,size(row,1))*m];
        end
        MarketMat = repmat(MarketMat,[J,1]);

        % Rankings post-counterfactual
        rankA_valns_sim = zeros(J,N/J,S);
        qualA_valns_sim = zeros(J,N/J,S);
        rankA_val_sim = zeros(J,N/J,S);
        qualA_val_sim = zeros(J,N/J,S);
        diffnsHOLD = zeros(N-1,S);
        diffHOLD = zeros(N-1,S);

        for s = 1:S
            shareIter = share0;
            shareKEEP = []; shareStKEEP = [];
            tol = 1; iter = 1;

            while tol > 1e-6 && iter < 300
                if cf == 11 || cf == 12 || cf == 13
                   ParcelIndex(:,cN) = Quality; 
                   ParcelUtility    = ParcelIndex*gamma';
                   FirmMatrix(:,end) = Quality; 
                   FirmUtility = FirmMatrix*delta';
                end
        
                % Sub-market profits and index
                ParcelUtility_sub = ParcelUtility + alpha*shareIter + ParcelError_Sim(1:N,s); 
                FirmUtility_sub   = FirmUtility   + FirmError_Sim(1:N,s); 

                Aoffer  = reshape(ParcelUtility_sub,[J,N/J]); 
                Baccept = reshape(FirmUtility_sub,[J,N/J]);

                % Sort preferences
                if out == 1
                    rankA = SortingMatInC_out(Aoffer); 
                    rankB = SortingResMatInC_out(Baccept);
                else
                    rankA = SortingMatInC(Aoffer); 
                    rankB = SortingResMatInC(Baccept);
                end

                % GS algorithm
                [~,FirmMatchOUT] = MatchAofferBaccept_GS_out(rankA,rankB,FirmMatchConstraint);
                [shareOUT,shareStOUT] = ShareMarket_st(FirmMatchOUT',MarketCount0,J);

                shareStKEEP = [shareStKEEP shareStOUT];
                shareKEEP   = [shareKEEP shareOUT];

                % tol = mean(abs(share0 - shareStOUT));
                if iter > 1
                %   tol = mean(abs(shareKEEP(:,iter - 1) - shareKEEP(:,iter)))
                    tol = sum(abs(shareKEEP(:,iter - 1) - shareKEEP(:,iter)));
                else
                %   tol = mean(abs(share0 - shareStOUT))
                    tol = sum(abs(shareIter - shareStOUT));
                end
                shareIter = shareStOUT;
            %       shareIter = mean(shareStKEEP,2);

                FirmMatchCF(1,1:size(FirmMatchOUT,2))  = FirmMatchOUT;
                shareOUTCF(1:size(shareOUT,1),1)       = shareOUT;
                iter = iter + 1;
            end
            FirmMatchCF(s,1:size(FirmMatchOUT,2))  = FirmMatchOUT;
            shareOUTCF_st(:,s)      = shareStOUT;
            QualityCF_st(:,s)       = Quality;

            [~,~,GpPMeanShare,GpFMeanShare,~,GpShareFirm,~,GpParcelCov,GpFirmCov] = GroupMoments_share2(ParcelIndex,FirmMatrix,FirmMatchOUT,J,MarketIndex0,MarketCount0,shareOUT);
            [~,~,PMean,FMean,JntDist,Cov,Pvar,Fvar,PFcovar,Pmean_firm,Pvar_firm,PFJntDist_firm,PFcovar_firm,~,~,~,~] = GroupMomentsSmall_new(ParcelIndex,FirmMatrix,FirmMatchOUT,J,MarketIndex0,MarketCount0,MomentSwitch);

            % Calculate the moments
            SampleMMean_cf(:,:,s)  = [FMean PMean];
            SampleMCov_cf(:,:,s)   = PFcovar;
            SampleMVar_cf(:,:,s)   = [Fvar Pvar];

            SampleJMean_cf(:,:,s)  = Pmean_firm;
            SampleJCov_cf(:,:,s)   = PFcovar_firm;
            SampleJVar_cf(:,:,s)   = Pvar_firm;
            SampleJDist_cf(:,:,s)  = PFJntDist_firm;

            SampleJntDist_cf(:,s)  = JntDist;
            SampleCov_cf(:,s)      = Cov;

            % Calculate share - Market
            Share_cf(:,s)          = shareOUT;
            ShareMean_cf(:,:,s)    = [GpFMeanShare GpPMeanShare];
            ShareFirm_cf(:,s)      = GpShareFirm;
            ShareCov_cf(:,:,s)     = [GpFirmCov GpParcelCov];

            % CF part II
            % Calculate the parcel and firm utility earned from each match across simulations
            ParcelIndex(:,cN) = QualityCF_st(:,s); 
            FirmMatrix(:,end) = QualityCF_st(:,s); 
            Quality_rsh       = reshape(QualityCF_st(:,s),[J,N/J]);

            Aoffer   = reshape(ParcelIndex*gamma' + alpha*shareOUTCF_st(1:N,s) + ParcelError_Sim(1:N,s),[J,N/J]);
            Baccept  = reshape(FirmMatrix*delta'  + FirmError_Sim(1:N,s),[J,N/J]);
            for i = 1:N/J
                if FirmMatchCF(s,i) > 0 
                    ParcelUtil_sim(s,i) = Aoffer(FirmMatchCF(s,i),i);
                    FirmUtil_sim(s,i)   = Baccept(FirmMatchCF(s,i),i);
                    Quality_sim(s,i)    = Quality_rsh(FirmMatchCF(s,i),i);
                    landsqft_sim(s,i)   = landsqft(FirmMatchCF(s,i),i);
                    gas_sim(s,i)        = gas(FirmMatchCF(s,i),i);
                    LandByMarket(MarketMat(1,i),1,s,v+1)   = LandByMarket(MarketMat(1,i),1,s,v+1) + landsqft(1,i);
                else
                    ParcelNotMatched_sim(s,i,1)   = mean(Baccept(:,i));
                    ParcelNotMatched_sim(s,i,2)   = std(Baccept(:,i));
                    ParcelNotMatched_sim(s,i,3)   = median(Baccept(:,i));
                    ParcelNotMatched_sim(s,i,4)   = min(Baccept(:,i));
                    ParcelNotMatched_sim(s,i,5)   = max(Baccept(:,i));

                    landsqft_notM_sim(s,i)      = landsqft(1,i);
                    gas_notM_sim(s,i)           = gas(1,i);
                end
                LandByMarket(MarketMat(1,i),2,s,v+1)  = LandByMarket(MarketMat(1,i),2,s,v+1) + landsqft(1,i);
            end
            
            % No value of concentration
            ParcelIndex(:,cN)  = QualityCF_st(:,s); 
            FirmMatrix(:,end)  = QualityCF_st(:,s); 
            Quality_rsh        = reshape(QualityCF_st(:,s),[J,N/J]);

            Aoffer   = reshape(ParcelIndex*gamma' - theta_score(1,cN - 1)*QualityCF_st(:,s) + ParcelError_Sim(1:N,s),[J,N/J]);
            rankA   = SortingMatInC(Aoffer); 
            for j = 1:J
                for i = 1:N/J
                    if rankA(j,i) > 0 
                    rankA_valns_sim(j,i,s) = Aoffer(j,rankA(j,i));
                    qualA_valns_sim(j,i,s) = Quality_rsh(j,rankA(j,i));
                    end
                end
            end

            % Add the effect of concentration
            Aoffer  = reshape(ParcelIndex*gamma' - theta_score(1,cN - 1)*QualityCF_st(:,s) + alpha*shareOUTCF_st(1:N,s) + ParcelError_Sim(1:N,s),[J,N/J]);
            rankA   = SortingMatInC(Aoffer);
            for j = 1:J
                for i = 1:N/J
                    if rankA(j,i) > 0 
                    rankA_val_sim(j,i,s) = Aoffer(j,rankA(j,i));
                    qualA_val_sim(j,i,s) = Quality_rsh(j,rankA(j,i));
                    end
                end
            end
           diff = rankA_valns_sim(:,:,s);
           diff = diff(:);
           diff = sort(diff,'descend');
           diff = [diff [diff(2:end); 0]];
           diffnsHOLD(:,s) = diff(1:end-1,1) - diff(1:end-1,2);
           diff = rankA_val_sim(:,:,s);
           diff = diff(:);
           diff = sort(diff,'descend');
           diff = [diff [diff(2:end); 0]];
           diffHOLD(:,s) = diff(1:end-1,1) - diff(1:end-1,2);             
        end
        
        shareOUTCF_st_mod(1:N,v+1)       = mean(shareOUTCF_st,2);
        QualityCF_mod(1:N,v+1)           = mean(QualityCF_st,2);
        
        SampleMMean_mod(1:M,1:JJ+NN,v+1) = mean(SampleMMean_cf,3);
        SampleMVar_mod(1:M,1:JJ+NN,v+1)  = mean(SampleMVar_cf,3);
        SampleMCov_mod(1:M,1:JJ*NN,v+1)  = mean(SampleMCov_cf,3);
        
        SampleJMean_mod(1:J,1:NN,v+1)    = mean(SampleJMean_cf,3);
        SampleJCov_mod(1:J,1:JJ*NN,v+1)  = mean(SampleJCov_cf,3);
        SampleJVar_mod(1:J,1:NN,v+1)     = mean(SampleJVar_cf,3);
        SampleJDist_mod(1:J,1:JJ*NN,v+1) = mean(SampleJDist_cf,3);
        
        SampleJntDist_mod(1:NN*JJ,v+1)   = mean(SampleJntDist_cf,2);
        SampleCov_mod(1:NN*JJ,v+1)       = mean(SampleCov_cf,2);

        Share_mod(1:M*J,v+1)             = mean(Share_cf,2);
        ShareMean_mod(1:M,1:NN+JJ,v+1)   = mean(ShareMean_cf,3);
        ShareFirm_mod(1:M,v+1)           = mean(ShareFirm_cf,2);
        ShareCov_mod(1:M,1:NN+JJ,v+1)    = mean(ShareCov_cf,3);

        ParcelUtil_mod(v+1,1:N/J)        = mean(ParcelUtil_sim,1);
        FirmUtil_mod(v+1,1:N/J)          = mean(FirmUtil_sim,1);
        Quality_mod(v+1,1:N/J)           = mean(Quality_sim,1);
        ParcelNotMatched_mod(1:N/J,:,v+1)= squeeze(mean(ParcelNotMatched_sim,1));

        landsqft_mod(v+1,1:N/J)          = mean(landsqft_sim,1);
        landsqft_notM_mod(v+1,1:N/J)     = mean(landsqft_notM_sim,1);
        gas_mod(v+1,1:N/J)               = mean(gas_sim,1);
        gas_notM_mod(v+1,1:N/J)          = mean(gas_notM_sim,1);

        rankA_valns_mod(1:J,1:N/J,v+1)   = mean(rankA_valns_sim,3);
        qualA_valns_mod(1:J,1:N/J,v+1)   = mean(qualA_valns_sim,3);
        rankA_val_mod(1:J,1:N/J,v+1)     = mean(rankA_val_sim,3);
        qualA_val_mod(1:J,1:N/J,v+1)     = mean(qualA_val_sim,3);
 
        diffnsMean(1:N-1,v+1)            = mean(diffnsHOLD,2);
        diffMean(1:N-1,v+1)              = mean(diffHOLD,2);

    end

end

clear Bads* app_val* prodTC* Legal* Extern* expFuture*
clear parcelTo* largeOp* prodwell* landman* Landman* landSizeTo* pipeDist* complaint* parcelInt*

clear data*

clear *_cf
clear PMean FMean JntDist Cov Pvar Fvar PFcovar Pmean_firm Pvar_firm PFJntDist_firm PFcovar_firm GpMean GpVar GpMean_byJ GpVar_byJ
clear GpPMeanShare GpFMeanShare GpShareMarket GpShareFirm GpFMean GpParcelCov GpFirmCov
clear FirmError_Sim ParcelError_Sim expFutureParcel* expFutureFirm*
clear Extern* prodTCwell* Legal* lateral* Landman* app_* env* complaint* Clause* Bonus* buffer* wellDist* pipe* prodwell* Bads*
clear shareStKEEP shareStOUT shareKEEP Aoffer Baccept
clear QualityCF_st shareOUTCF_st landSize* ToNorm*
clear shareIter share0
clear ParcelUtility* FirmUtility* diffHOLD diffnsHOLD data data0 rows row0 

save(['ModelSelect' num2str(datsize) '_sim' num2str(b) '_full_endqual_cf' num2str(cf) '.mat'])

